﻿using CsvHelper;
using Python.Runtime;
using System.Globalization;

namespace PythonnetDemo
{
    /*
    internal class Program1
    {
        const string PathToPythonDir = "D:\\Python311";
        const string DllOfPython = "python311.dll";

        static void Main(string[] args)
        {
            // 数据清洗
            // CleanData();
            // 傅里叶变换
            // FFT();
            // 低通滤波
            // ButterworthLowpass();
        }

        /// <summary>
        /// 傅里叶变换
        /// </summary>
        static void FFT()
        {

            try
            {
                Runtime.PythonDLL = Path.Combine(PathToPythonDir, DllOfPython);

                PythonEngine.Initialize();
                using (Py.GIL())
                {
                    dynamic pd = Py.Import("pandas");
                    dynamic np = Py.Import("numpy");
                    dynamic plt = Py.Import("matplotlib.pyplot");
                    dynamic fft = Py.Import("scipy.fftpack");

                    using dynamic scope = Py.CreateScope();
                    scope.Exec(@"def get_slice(net_array): return net_array[:len(net_array)//2]");

                    // 读取数据
                    var data = pd.read_csv("clean_data_row6.csv");
                    int listLength = data.__len__();
                    Console.WriteLine("读取长度:" + listLength);

                    // 进行傅里叶变换
                    var fft_result = fft.fft(data); // 对数据进行傅里叶变换
                    var frequencies = fft.fftfreq(listLength);

                    // 计算功率谱密度
                    var power_spectrum = np.square(np.abs(fft_result)) / listLength;
                    var frequencies_positive = scope.get_slice(frequencies);

                    double[] xAxis = frequencies.As<double[]>();
                    PyObject yAxisDatas = np.abs(fft_result);
                    double[][] yAxis = yAxisDatas.As<dynamic[]>()
                        .Select(s => (double[])s.As<double[]>())
                        .ToArray();

                    double[] xAxis2 = xAxis.Take(listLength / 2).ToArray();
                    PyObject yAxisDatas2 = power_spectrum;
                    double[][] yAxis2 = yAxisDatas2.As<dynamic[]>()
                        .Select(s => (double[])s.As<double[]>())
                        .Take(listLength / 2)
                        .ToArray();

                    ;

                    // 绘制频谱图和功率谱密度图
                    plt.figure(figsize: new dynamic[] { 12, 6 });

                    // 频谱图
                    plt.subplot(1, 2, 1);
                    plt.plot(frequencies, yAxisDatas);
                    plt.xlabel("Frequency (Hz)");
                    plt.ylabel("Amplitude");
                    plt.title("Frequency Spectrum");

                    // 功率谱密度图
                    plt.subplot(1, 2, 2);
                    plt.plot(frequencies_positive, scope.get_slice(power_spectrum));
                    plt.xlabel("Frequency (Hz)");
                    plt.ylabel("Power");
                    plt.title("Power Spectrum Density");

                    // 布局调整，防止重叠
                    plt.tight_layout();
                    // 显示图表
                    plt.show();
                }
            }
            catch (Exception e)
            {
                Console.WriteLine("报错了：" + e.Message + "\r\n" + e.StackTrace);
            }
        }

        /// <summary>
        /// Butter-worth低通滤波
        /// </summary>
        static void ButterworthLowpass()
        {
            int fs = 250;
            double filterCutoff = 50;
            int axis = 0;
            int factor = 0;
            int order = 8;

            try
            {
                // NOTE: set this based on your python install. this will resolve from
                // your PATH environment variable as well.
                Runtime.PythonDLL = Path.Combine(PathToPythonDir, DllOfPython);

                PythonEngine.Initialize();
                using (Py.GIL())
                {
                    dynamic pd = Py.Import("pandas");
                    dynamic np = Py.Import("numpy");
                    dynamic plt = Py.Import("matplotlib.pyplot");
                    dynamic signal = Py.Import("scipy.signal");

                    dynamic data = pd.read_csv("clean_data_row3.csv");
                    int listLength = data.__len__();
                    double wn = filterCutoff / (fs / 2.0);

                    PyTuple baTuple = signal.butter(order, wn, "low");
                    // PyTuple baTuple = signal.butter(4, 0.1, "low");
                    dynamic b = baTuple[0];
                    dynamic a = baTuple[1];
                    PyObject calcData = signal.filtfilt(b, a, data, axis: axis, padlen: listLength - 1);

                    double[][] calcDataArray = calcData.As<dynamic[]>().Select(s => (double[])s.As<double[]>()).ToArray();

                    plt.figure(figsize: new dynamic[] { 12, 6 });

                    // 频谱图
                    plt.subplot(1, 2, 1);
                    plt.plot(np.arange(listLength), calcData);
                    //plt.xlabel("");
                    //plt.ylabel("Amplitude");
                    plt.title("低通滤波");

                    // 布局调整，防止重叠
                    plt.tight_layout();
                    // 显示图表
                    plt.show();
                }
            }
            catch (Exception e)
            {
                Console.WriteLine("报错了：" + e.Message + "\r\n" + e.StackTrace);
            }
        }

        /// <summary>
        /// 数据清洗
        /// </summary>
        static void CleanData()
        {
            var originDatas = ReadCsvWithCsvHelper("data_row3.csv");

            var outliers = DetectOutliers(originDatas);

            var outlierHashset = new HashSet<int>(outliers);

            // 干净的数据
            var cleanDatas = originDatas.Where((r, index) => !outlierHashset.Contains(index)).ToList();

            try
            {
                // NOTE: set this based on your python install. this will resolve from
                // your PATH environment variable as well.
                Runtime.PythonDLL = Path.Combine(PathToPythonDir, DllOfPython);

                PythonEngine.Initialize();
                using (Py.GIL())
                {
                    dynamic pd = Py.Import("pandas");
                    dynamic np = Py.Import("numpy");
                    dynamic plt = Py.Import("matplotlib.pyplot");
                    dynamic fft = Py.Import("scipy.fftpack");

                    dynamic oData = np.array(originDatas.ToArray());
                    int oDataLength = oData.__len__();
                    dynamic data = np.array(cleanDatas.ToArray());
                    int dataLength = data.__len__();

                    // 绘制频谱图和功率谱密度图
                    plt.figure(figsize: new dynamic[] { 12, 6 });

                    // 频谱图
                    plt.subplot(1, 2, 1);
                    plt.plot(np.arange(oDataLength), oData);
                    plt.title("Original Datas");

                    // 功率谱密度图
                    plt.subplot(1, 2, 2);
                    plt.plot(np.arange(dataLength), data);
                    plt.title("Clean Datas");

                    // 布局调整，防止重叠
                    plt.tight_layout();
                    // 显示图表
                    plt.show();
                }
            }
            catch (Exception e)
            {
                Console.WriteLine("报错了：" + e.Message + "\r\n" + e.StackTrace);
            }
        }

        /// <summary>
        /// 检测异常值
        /// </summary>
        /// <param name="datas"></param>
        /// <returns>返回异常值在集合中的索引</returns>
        static List<int> DetectOutliers(List<double[]> datas)
        {
            List<int> outliers = new List<int>();
            var first = datas.First();
            for (int i = 0; i < first.Length; i++)
            {
                var values = datas.AsEnumerable().Select((row, index) => Tuple.Create(row[i], index)).ToArray();

                double q1 = Enumerable.OrderBy(values, x => x.Item1).ElementAt((int)(values.Length * 0.25)).Item1;
                double q3 = Enumerable.OrderBy(values, x => x.Item1).ElementAt((int)(values.Length * 0.75)).Item1;
                double iqr = q3 - q1;
                double lowerBound = q1 - 1.5 * iqr;
                double upperBound = q3 + 1.5 * iqr;

                outliers.AddRange(values.AsEnumerable()
                    .Where(row => row.Item1 < lowerBound || row.Item1 > upperBound)
                    .Select(row => row.Item2));
            }
            return outliers.Distinct().ToList();
        }

        /// <summary>
        /// 降采样，其实就是抽稀算法
        /// </summary>
        static List<double[]> DownSampling()
        {
            int factor = 2;
            int axis = 0;
            if (axis != 0 && axis != 1)
                throw new ArgumentException("Axis must be 0 or 1 for a 2D array.");

            var datas = ReadCsvWithCsvHelper("clean_data_row3.csv");

            int dim0 = datas.Count;
            var first = datas.First();
            int dim1 = first.Length;

            var result = new List<double[]>();
            if (axis == 0)
            {
                var xAxis = dim0 / factor;
                var yAxis = dim1;
                for (int i = 0; i < xAxis; i++)
                {
                    result.Add(datas[i * factor]);
                }
            }
            else if (axis == 1)
            {
                var xAxis = dim0;
                var yAxis = dim1 / factor;
                var item = new double[yAxis];
                for (int i = 0; i < xAxis; i++)
                {
                    var deviceData = datas[i];
                    for (int j = 0; j < yAxis; j++)
                    {
                        item[j] = deviceData[j * factor];
                    }
                    result.Add(item);
                }
            }
            return result;
        }

        /// <summary>
        /// 读取CSV数据
        /// </summary>
        /// <param name="filePath">文件路径</param>
        /// <returns></returns>
        static List<double[]> ReadCsvWithCsvHelper(string filePath)
        {
            
            // var result = new List<double[]>();
            // for (int i = 0; i < 10000; i++)
            // {
            //     result.Add([-0.01d]);
            // }
            // return result;
          
            using (var reader = new StreamReader(filePath))
            using (var csv = new CsvReader(reader, CultureInfo.InvariantCulture))
            {
                var result = new List<double[]>();
                // 如果你的CSV文件有标题行，可以调用ReadHeader来读取它们
                csv.Read();
                csv.ReadHeader();
                while (csv.Read())
                {
                    result.Add(new double[] {
                        csv.GetField<double>(0),
                        csv.GetField<double>(1),
                        csv.GetField<double>(2),
                        //csv.GetField<double>(3), 
                        //csv.GetField<double>(4), 
                        //csv.GetField<double>(5)
                    });
                }
                return result;
            }
        }
    }
    */
public class Program
{
    const string PathToPythonDir = "D:\\Python311";
    const string DllOfPython = "python311.dll";

    static void Main(string[] args)
    {
        // Butter-worth低通滤波
        ButterworthLowpass();
    }
    /// <summary>
    /// Butter-worth低通滤波
    /// </summary>
    static void ButterworthLowpass(int fs = 250, double filterCutoff = 50, int axis = 0, int order = 8)
    {
        try
        {
            Runtime.PythonDLL = Path.Combine(PathToPythonDir, DllOfPython);

            PythonEngine.Initialize();
            using (Py.GIL())
            {
                dynamic pd = Py.Import("pandas");
                dynamic np = Py.Import("numpy");
                dynamic plt = Py.Import("matplotlib.pyplot");
                dynamic signal = Py.Import("scipy.signal");

                dynamic data = pd.read_csv("clean_data_row3.csv");
                int listLength = data.__len__();
                double wn = filterCutoff / (fs / 2.0);

                PyTuple baTuple = signal.butter(order, wn, "low");
                dynamic b = baTuple[0];
                dynamic a = baTuple[1];
                PyObject calcData = signal.filtfilt(b, a, data, axis: axis, padlen: listLength - 1);
                int calcLength = data.__len__();
                double[][] calcDataArray = calcData.As<dynamic[]>().Select(s => (double[])s.As<double[]>()).ToArray();

                plt.figure(figsize: new dynamic[] { 12, 6 });

                // 原始数据
                plt.subplot(1, 2, 1);
                plt.plot(np.arange(listLength), data);
                plt.title("Original Datas");

                // 低通滤波
                plt.subplot(1, 2, 2);
                plt.plot(np.arange(calcLength), calcData);
                plt.title("Butterworth Lowpas Datas");

                // 布局调整，防止重叠
                plt.tight_layout();
                // 显示图表
                plt.show();
            }
        }
        catch (Exception e)
        {
            Console.WriteLine("报错了：" + e.Message + "\r\n" + e.StackTrace);
        }
    }

    /// <summary>
    /// 读取CSV数据
    /// </summary>
    /// <param name="filePath">文件路径</param>
    /// <returns>文件中数据集合，都是double类型</returns>
    static List<double[]> ReadCsvWithCsvHelper(string filePath)
    {
        using (var reader = new StreamReader(filePath))
        using (var csv = new CsvReader(reader, CultureInfo.InvariantCulture))
        {
            var result = new List<double[]>();
            // 如果你的CSV文件有标题行，可以调用ReadHeader来读取它们
            csv.Read();
            csv.ReadHeader();
            while (csv.Read())
            {
                result.Add(new double[] {
                    csv.GetField<double>(0),
                    csv.GetField<double>(1),
                    csv.GetField<double>(2),
                });
            }
            return result;
        }
    }
}
}

